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ABSTRACT 


The  superfluid  properties  of  disordered  double-layer  graphene  systems  are  in¬ 
vestigated  using  the  non-equilibrium  Green’s  function  (NEGF)  formalism.  The 
complexity  of  such  a  structure  makes  it  imperative  to  study  the  effects  of  lattice 
vacancies  which  will  inevitably  arise  during  fabrication.  Room-temperature  per¬ 
formance  characteristics  for  both  ideal  and  disordered  bilayer  graphene  systems 
are  compared  in  an  effort  to  illustrate  the  behavior  of  a  Bose-Einstein  condensate 
in  the  presence  of  lattice  defects  under  non-equilibrium  conditions.  The  study 
finds  that  lattice  vacancies  that  spread  throughout  the  top  layer  past  the  coherence 
length  have  a  reduced  effect  compared  to  the  ideal  case.  However,  vacancies  con¬ 
centrated  near  the  metal  contacts  within  the  coherence  length  significantly  alter 
the  interlayer  superfluid  transport  properties. 
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CHAPTER  1 


MOTIVATION 


From  high-performance  parallel  computing  to  hand-held  devices,  the  advent  of  ef¬ 
ficient  computer  processors  has  significantly  changed  the  face  of  the  technological 
landscape.  The  Internet  has  pervaded  almost  every  facet  of  business  information 
technology,  from  marketing  to  supply-chain  management.  None  of  this  revolution 
would  be  possible  without  mobile  and  efficient  data  processing.  Any  computer 
or  smart  phone  today  is  powered  by  complementary  metal- oxide- semiconductor 
(CMOS)  technology,  which  uses  gated  metal-oxide-semiconductor  field  effect 
transistors  (MOSFETs)  to  implement  Boolean  logic  gates  to  carry  out  any  par¬ 
ticular  computation.  For  the  better  part  of  40  years,  the  density  and  efficiency  of 
these  transistors  has  very  closely  followed  Moore’s  law,  which  states  the  number 
of  transistors  in  an  integrated  circuit  doubles  approximately  every  two  years  [1]. 
This  is  a  remarkable  trend  propelled  by  decades  of  innovation  in  fabrication  of 
silicon  CMOS  technology.  The  trend  is  expected  to  slow  in  subsequent  years, 
however,  as  the  transistor  sizes  and  densities  approach  fundamental  limits. 

Perhaps  the  most  insurmountable  of  these  limits  is  the  power  density  of  cur¬ 
rent  processors.  Heat  is  generated  when  current  flows  through  the  doped  silicon 
channel  of  a  MOSFET.  At  low  density,  the  heat  can  dissipate  with  simple  air  cool¬ 
ing  techniques.  Higher  transistor  densities,  however,  have  begun  to  heat  chips 
near  the  air  cooling  limit,  at  which  point  the  processor  would  become  too  hot  for 
normal  operation  and  fails.  At  the  same  time,  leakage  current  becomes  a  non- 
negligible  factor  in  undesirable  energy  consumption  as  transistors  shrink  to  di¬ 
mensions  where  the  effects  of  quantum  tunneling  through  the  gate  oxide  become 
apparent. 

A  large  variety  of  devices  have  been  proposed  to  overcome  these  inefficiencies. 
Three-dimensional  gate  structures  were  recently  introduced  in  fabrication  to  in¬ 
crease  channel  mobility  and  decrease  gate  voltages  in  silicon  channels.  Gallium 
arsenide  (GaAs)  and  other  III- V  channels  provide  an  increase  in  electron  mobility, 
decreasing  power  dissipation  and  the  resultant  heat  generation.  Organic  electron- 
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ics,  utilizing  graphene  or  carbon  nanotube  (CNT)  channels  have  more  recently 
been  studied  due  to  their  high  carrier  mobilities  and  thermal  conductivity.  Opto¬ 
electronic  devices,  using  photons  as  signal  carriers  rather  than  electrons,  are  also 
being  actively  pursued  and  have  found  practical  applications  in  large-scale  data 
transfer  utilizing  optical  fibers. 

While  this  research  has  provided  tremendous  insight  for  solutions  in  the  short 
run,  the  MOSFET’s  unavoidable  limitations  will  remain  a  problem  in  any  sys¬ 
tem  with  a  semiconducting  channel  due  to  dissipation  and  scattering  of  individual 
carriers  in  the  crystal  lattice.  A  superconducting  channel  would  resolve  these 
issues,  but  to  date  no  room-temperature  superconducting  system  has  been  real¬ 
ized,  and  the  channel  cannot  be  turned  on  or  off  via  a  gate  voltage.  Recent  work 
in  condensed  matter  physics,  however,  has  shown  evidence  of  other  condensed 
many-body  systems  that  can  flow  around  channel  defects  and  impurities  without 
dissipation  or  scattering  losses  in  the  same  fashion  as  a  superconductor.  In  Chapter 
2, 1  describe  this  many-body  phenomenon,  known  as  Bose-Einstein  condensation 
(BEC),  and  compare  it  to  conventional  superconductivity.  I  analyze  its  experi¬ 
mental  realization  in  a  variety  of  bosonic  systems  to  elucidate  how  condensation 
starkly  effects  various  macroscopic  behaviors.  In  Chapter  3, 1  further  discuss  the 
realization  of  superfluid  transport  in  coupled  quantum  well  systems,  and  extend 
this  idea  in  Chapter  4  to  explain  why  double-layer  graphene  systems  are  a  suitable 
candidate  for  room-temperature  operation  in  which  the  superfluid  can  be  switched 
on  or  off  with  a  double-gate  geometry.  In  Chapter  5, 1  provide  a  detailed  summary 
of  the  non-equilibrium  Green’s  function  (NEGF)  formalism  and  its  ability  to  de¬ 
termine  non-equilibrium  transport  dynamics  of  a  system.  In  Chapter  6,  I  utilize 
NEGF  to  predict  the  observable  behavior  of  the  double-layer  graphene  system  in 
both  ideal  and  disordered  systems  in  order  to  understand  how  imperfect  graphene 
channels  degrade  transport  dynamics. 
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CHAPTER  2 


INTRODUCTION 


2.1  Superconductivity 

A  boson  is  defined  as  any  quasiparticle  with  integer  spin.  Bose  first  showed  in 
1925  that  photons,  spin-0  light  quanta,  adhere  to  statistics  qualitatively  different 
than  fermions,  particles  with  half-integer  spin.  Einstein  subsequently  extended 
Bose’s  study  to  matter  to  show  that  all  integer-spin  particles  obey  Bose-Einstein 
statistics,  where  the  expected  number  of  particles  in  a  non-degenerate  energy  state 
i  is 

Ui  =  e(Ei-n)/kBT  _  X 

where  E*  is  the  energy  of  state  i,  /i  is  the  chemical  potential,  T  is  the  temperature 
and  ks  is  Boltzmann’s  constant.  At  energies  well  below  the  chemical  potential, 
the  expected  number  of  particles  becomes  much  larger  than  1.  This  is  in  stark  con¬ 
trast  to  Fermi-Dirac  statistics  for  fermions,  or  spin-|  particles  like  electrons  and 
holes.  The  Pauli  exclusion  principle  prevents  fermions  from  occupying  the  same 
quantum  state,  limiting  the  expectation  number  to  less  than  1  for  every  fermion 
quantum  state. 

At  low  temperature,  a  large  fraction  of  the  bosons  in  a  system  decay  to  the 
lowest  quantum  state.  Condensation  in  the  lowest  quantum  state  manifests  itself 
with  a  variety  of  phenomena  on  the  macroscopic  scale,  including  superconductiv¬ 
ity.  The  Bardeen-Cooper-Schrieffer  (BCS)  theory  of  superconductivity  [2]  was  a 
crowning  achievement  of  quantum  mechanics  in  condensed  matter  systems.  The 
discovery,  now  over  fifty  years  old,  explains  superconductivity  as  the  result  of  con¬ 
densation  of  composite  bosons  into  a  ground  state.  BCS  predicted,  however,  that 
two  electrons  can  couple  to  form  a  composite  boson  quasiparticle  called  a  Cooper 
pair.  Electrons  are  repelled  from  one  another  by  the  Coulomb  interaction,  but  at 
very  low  temperature  a  phonon-mediated  electronic  attraction  can  overcome  the 
Coulomb  repulsion.  The  resulting  Cooper  pair  has  an  energy  lower  than  the  en- 
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ergy  of  a  bare  electron,  and  opens  a  gap  about  the  Fermi  energy,  EF.  Cooper  pairs 
are  capable  of  flowing  through  the  superconducting  channel  without  dissipation, 
a  macroscopic  realization  of  the  many-body  coherence.  The  pairing  is  very  weak, 
however,  and  superconductivity  is  quickly  lost  as  temperature  increases  above  the 
critical  temperature  Tc  -  the  temperature  at  which  phonons,  vibrations  in  the  crys¬ 
tal  lattice,  dominate  electron  motion  and  annihilate  the  weak  Cooper  pairs. 

While  the  discovery  of  superconductivity  has  bred  a  wealth  of  research  into 
high  Tc  unconventional  superconductors,  whose  critical  temperatures  approach 
150  K  [3],  Tc  remains  well  below  room  temperature  due  to  the  long-range  phonon- 
mediated  nature  of  the  condensation.  However,  if  a  strong  Coulomb  attraction 
between  oppositely  charged  particles  replaced  the  weak  attraction  in  supercon¬ 
ductivity,  the  bosons  could  condense  and  survive  at  much  higher  temperatures.  In 
the  following  section,  I  will  describe  such  an  interaction  and  show  how  the  pairing 
can  be  large  enough  to  remain  condensed  at  room  temperature. 


2.2  Bose-Einstein  Condensation 


Although  BEC  was  first  predicted  in  the  1920s,  experimental  realization  did  not 
occur  until  1995  with  a  gas  of  weakly  interacting  rubidium  atoms  cooled  to  nearly 
20  nK  [4].  Anderson  et  al.  found  a  signature  of  BEC  when  the  concentration  of 
rubidium  atoms  in  the  zero-velocity  ground  state  increased  abruptly  as  the  sample 
was  lowered  below  the  critical  temperature,  which  for  a  three-dimensional  non¬ 
interacting  gas  is 


T 

±  C 


2nh2  f  n 
rnRbkB  V C (3/2) 


2/3 

«  50  nK. 


(2.2) 


where  mRb  is  the  mass  of  a  rubidium  atom,  n  is  the  number  of  particles  in  the 
system,  and  (  is  the  Riemann  zeta  function. 

While  the  rubidium  atom  experiment  provided  definitive  evidence  of  BEC,  a 
sub-micro  Kelvin  critical  temperature  is  of  little  practical  value.  However,  this 
realization  is  in  the  weakly  interacting  regime  with  very  heavy  atoms  (note  that 
critical  temperature  is  inversely  proportional  to  mass  in  Eq.  2.2).  Subsequent  re¬ 
search  has  shown  that  condensates  can  be  formed  with  weakly  interacting  exciton- 
polariton  quasiparticles,  which  are  created  with  a  continuous  wave  photo-pump  in 
a  semiconductor  cavity.  As  long  as  the  polariton  generation  rate  is  larger  than  the 
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rate  at  which  the  quasiparticles  recombine,  a  condensate  can  persist.  Exciton- 
polariton  quasiparticles  have  effective  masses  109  times  lighter  than  rubidium 
atoms,  so  that  they  are  suitable  for  observed  condensation  at  19  K  [5;  6].  Further 
research  has  shown  that  polaritons  with  even  smaller  effective  masses  can  have 
critical  temperatures  above  room  temperature  in  pristine  bulk  GaN  microcavities 
[7;  8], 

Polariton  condensates  provide  a  unique  environment  to  observe  various  signa¬ 
tures  of  BEC.  The  most  important  of  these  for  the  scope  of  this  thesis  is  dissipa¬ 
tionless  flow,  or  zero  viscosity,  known  as  superfluidity.  The  BEC  is  a  superfluid 
state,  meaning  it  can  flow  freely  and  without  diffusion  through  a  system  in  the 
same  way  a  Cooper  pair  flows  freely  through  a  superconducting  channel.  The 
best  known  observation  of  this  is  the  ability  of  superfluid  He4  to  escape  out  of  a 
container  by  flowing  along  the  walls  of  the  container  [9].  In  a  similar  fashion, 
Amo  et  al.  [10]  showed  that  droplets  of  polariton  condensation  can  be  excited  to 
flow  through  the  channel  in  which  they  are  formed  with  a  simple  light  pulse.  The 
droplets  move  with  a  constant  high  group  velocity,  roughly  1.0  x  106  m/s,  and 
are  able  to  move  around  defects  without  scattering.  When  the  polariton  superfluid 
reaches  a  defect  with  size  comparable  to  its  own,  it  tends  to  split  the  superfluid 
into  two  droplets,  but  phase  coherence  is  not  lost  and  group  velocity  remains  un¬ 
perturbed.  Although  the  exciton-polariton  is  charge  neutral,  the  observed  trans¬ 
port  dynamics  offer  strong  insight  to  the  technological  ramifications  of  superfluid 
transport  -  macroscopic  coherence  that  precludes  diffusion  from  material  defects. 
In  the  next  section,  I  will  describe  a  system  in  which  charge  can  be  separately 
transfered  in  a  excitonic  superfluid  system,  as  well  as  experimental  realization  in 
the  quantum  Hall  regime. 
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CHAPTER  3 


TRANSPORT  IN  EXCITONIC 
SUPERFLUIDS 


3.1  Crossed  Andreev  Reflection 

From  a  theoretical  standpoint,  an  indirectly  bound  exciton  superfluid  can  form 
when  two  channels  are  stacked  in  close  proximity  to  one  another  and  doped  so  that 
a  large  number  of  electrons  (holes)  occupy  the  top  (bottom)  layer.  The  carriers  in 
the  top  layer  will  form  indirect  excitons  (bosons)  with  holes  in  the  bottom  layer 
via  an  interlayer  Coulomb  interaction,  but  are  prevented  from  recombining  by  a 
small  spacer  dielectric.  This  is  analogous  to  the  polariton  condensates  discussed  in 
the  previous  section,  with  the  exception  that  photo-pumping  is  not  necessary  for 
condensation,  and  the  constituent  components  of  the  condensate  are  separately 
contacted. 

It  is  vital  to  first  understand  the  transport  properties  that  can  be  achieved  in  such 
an  indirectly  bound  excitonic  superfluid  system.  To  obtain  a  more  in-depth  under¬ 
standing  of  the  interlayer  transport  in  double-layer  graphene  systems,  I  first  draw 
a  comparison  to  the  Andreev  reflections  [12]  that  occur  at  a  metal- superconductor 
interface.  When  an  electron  with  energy  less  than  the  superconducting  band  gap  is 
injected  into  a  superconducting  system,  the  injected  electron  penetrates  a  certain 
distance  before  propelling  a  Cooper  pair  across  the  superconducting  channel,  as 
seen  in  Fig.  3.1a.  A  hole  of  opposite  spin  is  retro-reflected  back  into  the  metal  con¬ 
tact  to  conserve  momentum.  Similar  physics  occur  when  an  electron  with  energy 
less  than  the  superfluid  band  gap  is  injected  into  the  excitonic  superfluid  channel. 
Figure  3.1b  shows  that  the  injected  electron  will  penetrate  a  short  distance  in  the 
top  layer  before  driving  an  exciton  that  moves  across  the  channel.  This  process 
corresponds  to  an  electron  being  injected  into  the  top  left  contact  that  causes  an 
exciton  to  move  across  the  system  while  an  electron  is  reflected  into  the  bottom 
left  contact  to  conserve  current  [11].  The  injected  and  retro-reflected  currents 
are  separated  in  the  top  left  and  bottom  left  contacts,  and  the  resulting  physics  is 
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(b) 

Normal  Superfluid  Normal 


Figure  3.1:  (a)  Depiction  of  Andreev  reflection  in  metal- superconductor- metal 
system.  Non-equilibrium  contact  bias  will  inject  an  electron  of  a  certain  spin 
(black  dot)  into  the  system  from  the  left  contact.  The  action  will  cause  a  Cooper 
pair  to  flow  along  the  channel,  resulting  in  a  single  electron  coming  out  of  the 
right  contact  and  a  hole  (red  dot)  of  opposite  spin  from  the  left  contact  to  conserve 
momentum,  (b)  Depiction  of  crossed  Andreev  reflection  in  an  indirectly  bound 
exciton  superfluid.  When  an  electron  is  injected  in  the  top  left  contact  in  non¬ 
equilibrium,  an  exciton  will  be  driven  along  the  channel,  resulting  in  an  electron 
coming  out  of  the  top  right  contact  and  bottom  left  contact  to  conserve  momentum. 
Image  from  [11]. 
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termed  crossed  Andreev  reflection. 

The  drag-counterflow  geometry  (Vtl  =  — Vtr ;  Vbl  =  Vrr  =  0  V)  will  cause 
an  electron  to  be  reflected  into  the  bottom  left  contact  and  a  hole  into  the  bottom 
right  contact  [13],  inducing  an  effective  current  through  the  bottom  layer  only 
when  a  condensate  exists.  All  non-equilibrium  configurations  described  below 
are  in  this  geometry.  The  entire  process  results  in  a  condensate  current,  due  to  the 
propagating  exciton,  and  a  quasiparticle  current,  caused  by  the  injected  and  retro- 
reflected  individual  carriers.  The  quasiparticle  current  is  only  nonzero  within  the 
coherence  length  (Lc),  which  is  the  maximum  length  an  injected  particle  pene¬ 
trates  into  the  superfluid  gap  before  triggering  the  exciton  scattering  event.  While 
the  microscopic  attributes  of  the  superfluid  and  the  dynamics  of  the  Andreev  re¬ 
flection  may  not  be  readily  apparent  in  experiments,  they  will  have  a  macroscopic 
effect  on  the  observable  interlayer  current.  The  next  section  will  discuss  the  real¬ 
ization  of  such  a  system  using  coupled  quantum  wells  in  the  quantum  Hall  regime. 


3.2  Experimental  Evidence 

The  quantum  Hall  effect  (QHE),  first  shown  experimentally  in  1980  [14],  was 
a  landmark  discovery  in  condensed  matter  physics  for  its  ability  to  determine  the 
fine- structure  constant,  a,  and  the  quantum  of  conductance  with  high  precision.  In 
the  QHE,  a  strong  magnetic  field  of  15  T  is  applied  perpendicularly  to  a  pristine 
semiconductor  quantum  well  which  can  be  treated  as  a  2-dimensional  electron 
gas  (2DEG).  At  extremely  low  temperatures,  the  magnetic  field  dominates  the 
electron  dynamics  so  that  electrons  occupy  quantized  cyclotron  orbits  known  as 
Landau  levels  with  cyclotron  frequency  uic  =  qB /m,  where  B  is  magnetic  field 
strength,  q  is  the  charge  of  the  electron,  and  m  is  its  effective  mass.  Electrons  in 
each  Landau  level  have  discrete  energy  values, 

1 

En  =  huc(n  +  -),  (3.1) 

which  are  analogous  to  the  quantum  harmonic  oscillator.  At  weak  magnetic  fields, 
Landau  quantization  is  responsible  for  Shubnikov-de  Haas  oscillations,  a  wll- 
known  relation  in  which  the  longitudinal  resistance  oscillates  as  the  magnetic 
field  is  modulated.  Stronger  magnetic  fields  further  separate  the  discrete  Landau 
energies,  and  as  more  Landau  levels  are  pushed  above  Ep,  the  Hall  conductivity 
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decreases  according  to 

q 2 

°xy  =  z'-p  (3.2) 

where  v  is  the  number  of  Landau  levels  below  Ep.  This  is  due  to  the  fact  that  only 
edge  states  are  available  for  transport  in  the  quantum  Hall  regime;  longitudinal 
transport  is  prevented  by  the  confining  magnetic  field  unless  a  Landau  level  is 
crossing  E^.  The  QHE  thus  shows  that  conductance,  like  so  many  other  physical 
observables,  takes  on  a  quantized  value  proportional  to  some  integer  times  the 
quantum  of  conductance,  q2/h,  and  the  number  of  filled  Landau  levels  can  be 
extracted  easily  from  the  conductance.  It  is  important  to  note  that  conductance 
does  not  depend  on  semiconductor  material  or  length  because  carrier  edge  states 
are  completely  protected  from  backscattering  in  the  quantum  Hall  regime. 

In  recent  years,  significant  experimental  [15-18]  and  theoretical  [19-23]  progress 
has  been  made  in  the  search  for  BEC  in  coupled  quantum  wells  in  the  quantum 
Hall  regime.  Consider  two  GaAs  quantum  wells  separated  by  a  thin  insulating  bar¬ 
rier  layer.  When  the  individual  filling  factor  of  each  2DEG  is  tuned  to  v  =  1/2, 
so  that  the  total  filling  factor  is  v total  =  1,  the  density  of  electrons  and  holes  in 
the  lowest  Landau  level  is  equivalent.  The  system  thus  develops  a  spontaneous 
interlayer  phase  coherence  when  the  separation  between  the  two  layers  is  small 
enough  to  allow  interlayer  Coulomb  pairing,  which  occurs  at  less  than  roughly  10 
nm  [24].  The  holes  (electrons)  in  one  layer  form  indirect  excitons  with  the  elec¬ 
trons  (holes)  in  the  opposite  layer.  The  layer  separation  greatly  increases  exciton 
lifetimes,  which  otherwise  would  be  too  small  for  conventional  excitonic  conden¬ 
sation  in  bulk  semiconductors  [25]  because  the  electron  and  hole  are  not  hindered 
from  immediately  recombining.  Thus  the  layer  separation  is  finely  tuned  so  that 
it  is  small  enough  for  Coulomb  attraction  and  exciton  formation  but  large  enough 
to  prevent  tunneling  and  recombination. 

The  indirect  excitons,  being  composite  bosons,  form  a  BEC  in  the  same  fashion 
as  the  rubidium  atoms  and  exciton-polaritons.  When  carriers  are  injected  into  one 
of  the  two  2DEGs,  excitons  in  the  condensate  are  propagated  along  the  channel, 
the  physics  of  which  is  known  as  Andreev  reflection  and  is  discussed  at  length  in 
Section  3.1.  Thus  a  current  can  be  observed  in  the  opposite  layer  in  spite  of  a  bot¬ 
tom  layer  contact  bias.  This  is  known  as  resonant  tunneling  current,  which  peaks 
at  a  filling  factor  of  exactly  one-half  [26] .  Tunneling  current  is  linearly  propor¬ 
tional  to  the  total  current  in  the  system.  Critical  tunneling  current,  the  maximum 
attainable  current  in  the  BEC  phase,  was  observed  to  be  1.5  nA  [27]  before  the 
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condensate  is  destroyed  by  the  influx  of  individual  quasiparticles,  at  which  point 
tunneling  becomes  negligible. 

Critical  current  is  extremely  sensitive  to  filling  factor,  and  quickly  drops  to  a 
negligible  amount  when  total  filling  factor  is  below  0.9  or  above  1.1  [28].  This 
lends  credence  to  the  theory  that  a  BEC  can  only  form  in  a  coupled  quantum 
well  system  when  electron  and  hole  quasiparticle  densities  are  equivalent  in  each 
layer.  Inequalities  in  the  quasiparticle  concentrations  hinder  exciton  formation  by 
screening  out  interlayer  interaction.  The  resonant  tunneling  phenomenon  is  lost  at 
temperatures  outside  of  the  quantum  Hall  regime  [26].  Although  the  critical  tem¬ 
perature  for  excitonic  condensation  is  high,  equal  carrier  concentrations  can  only 
be  achieved  in  the  quantum  Hall  regime,  which  requires  sub-Kelvin  temperatures 
and  magnetic  fields  on  the  order  of  10  T.  Recent  theoretical  studies,  however,  have 
predicted  this  can  be  resolved  by  replacing  the  coupled  semiconductor  quantum 
well  heterostructure  with  a  double-layer  graphene  system.  In  the  next  section,  I 
will  analyze  carrier  transport  in  graphene  and  show  how  its  properties  make  it  an 
ideal  candidate  for  indirectly  bound  excitonic  condensation  with  critical  tempera¬ 
tures  approaching  room  temperature. 
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CHAPTER  4 


GRAPHENE 


4. 1  Electronic  Properties  of  Graphene 

The  unique  electronic  properties  of  graphene  -  an  individual  sheet  of  graphite  - 
were  first  predicted  in  1946,  long  before  its  experimental  discovery  [29].  Figure 
4. 1  shows  a  top-down  view  of  the  hexagonal  arrangement  of  the  carbon  atoms  in 
real  space  in  the  graphene  sheet  as  well  as  the  first  Brillouin  zone  in  reciprocal 
space.  The  basis  vectors  are 


a  i  —  I a0 


0*2  —  \a0\ 


(4.1) 


where  the  magnitude  of  the  basis  vectors  is  |a0|  =  V3acc.  The  distance  between 
nearest  neighbor  carbon  atoms  is  acc  rD.142  nm. 

In  order  to  calculate  the  band  structure  of  graphene,  a  tight  binding  model  is 
used  in  which  a  single  pz  orbital  per  carbon  atom  interacts  with  nearest  neighbor 
sites.  This  approximation  works  well  at  low  energies  near  the  Fermi  energy;  which 
is  the  only  important  energy  range  for  electrical  conduction,  because  other  orbitals 
are  either  well  above  or  well  below  Ep.  The  2x2  Hamiltonian  can  thus  be  written 
in  reciprocal  space  by  summing  over  any  unit  cell  and  its  six  neighboring  unit 
cells, 


H(k)  =  -t 


Q  j  gik-ai  j  eik-a2  | 

■j^  j  g  ik-ai  j  g  ik‘Q,2  j  g  ik'CLQ  q 


(4.2) 

where  r=— 3.03  eV  is  the  bonding  energy  between  carbon  atoms  and  (2  3  —  a  i  —  &2  — 
a0y.  The  dispersion  relation,  E (k),  can  be  calculated  by  solving  for  the  eigenvalues 
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(a) 


(b) 


Figure  4.1:  (a)  Real  space  depiction  of  a  graphene  lattice.  The  unit  cell  encom¬ 
passes  two  carbon  atoms,  with  basis  vectors  a\  and  a2.  (b)  The  first  Brillouin 
zone  of  graphene  in  reciprocal  space,  with  basis  vectors  b\  and  b2.  The  K  and  K’ 
points  in  the  Brillouin  zone  are  responsible  for  massless  transport  at  low  energy. 
Adapted  from  [30]. 


of  Eq.  4.2, 

E{k)  =  ±.t\J 3  +  2  cos  k  ■  ai  +  2  cos  k  ■  a2  +  2  cos  k  ■  a^,  (4.3) 

where  positive  (negative)  energies  denote  the  conduction  (valence)  bands  for  the 
graphene  lattice.  Graphene  is  thus  particle-hole  symmetric  about  Ep,  an  important 
property  which  will  be  exploited  in  Section  4.2  to  get  the  equivalent  quasiparticle 
densities  necessary  for  exciton  condensation.  This  property  is  in  stark  contrast  to 
silicon  and  other  bulk  semiconductors,  which  have  conduction  band  minima  and 
valence  band  maxima  that  do  not  lie  at  the  same  energy  or  wave  vector  and  are 
not  particle-hole  symmetric. 
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The  basis  vectors  in  reciprocal  space  must,  for  any  system,  satisfy 


*2716 ij , 


(4.4) 


where  StJ  is  the  Kronecker  delta.  For  graphene,  the  reciprocal  space  basis  vectors 
are 


4n 

ci  o\/3 


b2 


4tt  A  \[?y  \ 

AT!  \2X  ~  ~YyJ 


(4.5) 


to  satisfy  Equations  4.1  and  4.4.  The  K  and  K’  valleys  of  the  Brillouin  zone, 
shown  in  Fig.  4.1b,  are  located  at  wave  vector 


(n±^)bi  +  (m  +  ^2), 


n,  m  £  7L 


(4.6) 


where  n  and  m  are  integers.  When  k k  is  plugged  into  the  dispersion  relation 
equation,  the  energy  vanishes, 


2  2  4 

E(kx)  =  ±t\/  3  +  2  cos  ±-7r  +  2  cos  +  -7T  +  2  cos  ±-7r  =  0.  (4.7) 

o  o  o 


When  solved  numerically,  E (k)  follows  a  linear  dispersion  about  each  K/K’  point 
of  the  Brillouin  zone,  as  seen  in  Fig.  4.2.  The  linear  dispersion  is  symmetric 
about  the  Dirac  point,  ED,  in  the  approximate  energy  range  of  ED  ±  1  eV.  The 
Fermi  velocity  is  constant  in  the  linear  regime  and  proportional  to  the  slope  of  the 
dispersion, 

1  dE 

vF=^—tt  106m/s,  (4.8) 

which  is  roughly  1/300  the  speed  of  light  in  a  vacuum.  The  linear  band  structure 
means  that  carriers  act  as  massless  Dirac  fermions  in  graphene,  because  effective 
mass  is  proportional  to  d2E/dk2. 

The  electrical  properties  of  graphene  are  thus  unique  and  quite  remarkable.  Un¬ 
like  as  in  any  semiconducting  material,  electron  and  hole  quasiparticles  are  able 
to  flow  in  bulk  graphene  at  equivalent  velocities  due  to  their  equivalent  mobilities. 
For  this  reason,  graphene  has  become  an  area  of  immense  study  since  its  initial 
experimental  discovery  [31].  Proof  of  the  high  carrier  mobility  was  seen  from 
observation  of  the  QHE  in  graphene  at  room  temperature  [32].  The  predicted 
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Figure  4.2:  The  dispersion  relation,  E(k),  of  graphene  is  isotropic  and  linear  at 
low  energy.  Zero  energy  modes  exist  at  the  K  and  K’  points  of  the  Brillouin  zone. 
There  are  a  total  of  six  such  points  in  the  first  Brillouin  zone,  but  only  two  are 
independent.  Image  from  [33]. 


Fermi  velocity,  is  extremely  high,  and  allows  pristine  graphene  samples  to 
have  equivalent  room  temperature  carrier  mobilities  of  2  x  105  cnrV  's_l,  sev¬ 
eral  orders  higher  than  bulk  silicon  (<1400  cm2V_1s_1  [34]).  In  Section  4.2,  it  is 
shown  that  these  properties  make  graphene  an  ideal  choice  for  double-layer  exci- 
tonic  superfluidity.  Unfortunately,  efficient  graphene  growth  has  proven  difficult, 
and  the  best  large-scale  synthesis  techniques,  performed  by  epitaxial  growth  on 
copper  foils,  yield  carrier  mobilities  several  orders  smaller  than  the  best  samples 
[35].  As  such,  disorder  in  the  graphene  channel  becomes  an  important  factor  in 
superfluid  transport,  and  is  discussed  in  the  superfluid  regime  in  Section  6.3. 


4.2  Indirect  Condensation  in  Double-Layer  Graphene 

In  Chapter  3,  it  is  shown  that  a  BEC  can  arise  in  coupled  quantum  well  systems 
in  the  quantum  Hall  regime  when  each  layer  has  a  filling  factor  of  v  =  |.  The 
quantum  Hall  regime  is  necessary  for  particle-hole  symmetry  about  the  Fermi  en¬ 
ergy,  so  that  each  layer  has  an  equivalent  quasiparticle  density.  The  system  is  not 
suited  for  room-temperature  application  because  the  QHE  requires  very  large  and 
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precise  magnetic  fields  and  very  low  temperatures.  However,  recent  theoretical 
studies  [36;  37]  have  predicted  the  the  same  effect  is  possible  in  spatially  sepa¬ 
rated  monolayers  of  graphene,  as  depicted  in  Fig.  4.3.  Mean-field  calculations 
have  shown  that,  when  the  two  monolayers  are  separated  by  a  Si02  spacer  dielec¬ 
tric  on  the  order  of  1  nm  thick  and  the  Fermi  energy  of  the  top  (bottom)  monolayer 
is  gated  to  ED  +  0.4  eV  (ED  -  0.4  eV),  Tc  can  approach  room  temperature  [38]. 
In  the  next  chapter,  I  will  explain  the  NEGF  formalism  which  is  implemented  to 
study  the  transport  dynamics  of  the  double-layer  graphene  system. 


Figure  4.3:  Schematic  depiction  of  the  double-layer  graphene  system.  A  1  nm 
SiOo  dielectric  separates  the  two  graphene  monolayers,  whose  Fermi  levels  are 
set  to  ±0.4  eV  by  the  top  and  bottom  gates,  Vtg  and  V BG,  to  effect  equivalent  but 
oppositely  charged  quasiparticle  densities  in  each  layer.  Four  contacts  on  either 
side  of  the  monolayers  are  responsible  for  inducing  a  counterflow  current. 
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CHAPTER  5 


NON-EQUILIBRIUM  GREEN’S 
FUNCTION  FORMALISM 


5.1  Introduction 

Figure  5.1  shows  a  generic  system  to  be  simulated  with  NEGF  [39;  40].  The 
electron  dynamics  of  the  channel  are  completely  described  by  its  Hamiltonian,  H, 
which  is  connected  to  two  contact  reservoirs  with  respective  chemical  potential 
/r L  and  nR.  Contacts  couple  to  the  channel  via  self-energy  matrices,  £ L  and  £/,>, 
which  will  be  defined  in  Section  5.2,  to  induce  a  non-equilibrium  current  flow. 
The  simulation  assumes  ballistic  transport,  so  self-energy  from  scattering,  £5,  is 
neglected.  The  system  is  insulated  from  gate  contacts,  whose  effects  are  realized 
via  a  three-dimensional  Poisson  solver,  which  solves  for  the  potential  at  a  given 
point  in  real  space  based  on  the  electron  density, 

V2<Hf)  =  ^[ND(f)  -  n(r)],  (5.1) 

where  o{r)  is  the  potential  profile  and  ND ( r)  is  the  carrier  density  from  donors, 
which  is  nonzero  only  for  doped  semiconductors.  The  electron  density,  n(r),  is 
evaluated  in  the  NEGF  kernel  (see  Eq.  5.15),  thus  the  NEGF  and  Poisson  algo¬ 
rithms  must  be  repeated  until  a  self-consistent  solution  is  achieved.  The  system 
is  easily  expandable  to  systems  with  several  contacts  or  gate  configurations  by 
adding  additional  self-energy  matrices  to  NEGF  or  additional  sources  to  the  Pois¬ 
son  solver.  In  the  case  of  Fig.  4.3,  there  are  four  contacts  connected  to  the  system 
and  two  gates  electrostatically  dope  the  monolayers. 

The  channel  is  defined  to  be  three-dimensional  in  nature,  with  Nr,  Ny,  and  Nz 
lattice  points  in  the  x,  y,  and  z  directions.  The  full  three-dimensional  real  space 
Hamiltonian  matrix  can  be  written  in  terms  of  its  directional  hopping  amplitudes 
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Figure  5.1:  A  generalized  NEGF  system.  Two  contacts  with  chemical  potentials 
fiL  and  ///,.  place  the  system  with  Hamiltonian  H  in  a  non-equilibrium  state  to 
induce  current  flow.  An  assortment  of  gates  with  potential  V(;  electrostatically 
dopes  the  system  according  to  the  Poisson  equation. 


in  the  following  nested  block  tridiagonal  structure, 


H0 

Hx 

0 

0 

0 

H| 

Ho 

H, 

0 

0 

0 

HI 

H0 

H, 

0 

0 

0 

0 

0 

0 

Hi 

H0 

H, 

0 

0 

0 

HI 

Ho 

(5.2) 


In  Equation  5.2,  Hi  =  [I(atxath)  <8>  Hopz\  where  \nxnv)  's  the  identity  matrix 
of  dimension  NrNy,  Hop^  is  the  hopping  amplitude  in  the  z  direction,  and  <8) 
is  the  Kronecker  product.  The  on-diagonal  component  is  itself  a  nested  block 
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tridiagonal  matrix, 


o 

o 

X 
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0 

Ho, 

o 

o 

X 

(5.3) 


The  nested  matrix  H0  similarly  has  components  H0i  =  ll(Ny)  <S)  Hopx ]  and 


AA 

HopY 

0 

0 

0 

Ho'Py 

AA 

Hopy 

0 

0 

0 

Hop\r 

AA 

Hopy 

0 

0 

0 

0 

0 

0 

H0Py 

AA 

Hopy 

0 

0 

0 

Hop'y 

AA 

(5.4) 


Note  that  H  and  H0  are  block  tridiagonal,  but  H0o  is  not  when  periodic  boundary 
conditions  in  y  are  used. 


5.2  Contact  Self-Energies 

The  channel  is  only  connected  to  the  contacts  along  the  width  at  x  =  1  and  x  = 
Nx,  for  various  points  in  z,  so  that  the  self-energy  matrices  of  the  left  and  right 
contacts,  5^  and  S2,  have  a  nonzero  component  dimension  of  Ny  x  Ny.  The  full 
self-energy  matrix  has  dimensions  equivalent  to  the  full  Hamiltonian,  but  is  zero  at 
all  points  not  connected  to  the  lead.  The  contacts  are  assumed  to  be  semi-infinite 
so  that  the  surface  Green’s  function,  g)  ,  self-consistently  satisfies  the  relation 

g 8l(E)  =  [(E  +  iO+)I(JVl/)  -  H00  +  H^g£H0i].  (5.5) 
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The  nonzero  component  of  the  self-energy  matrix  can  then  be  computed  as 

Sr  (E)  =  t4g£(£)H01.  (5.6) 

£2(E)  can  be  calculated  in  the  same  fashion;  the  nonzero  components  of  the  self¬ 
energies  are  equivalent  in  this  particular  case.  The  imaginary  part  of  the  self¬ 
energy  matrix  is  the  energy  level  broadening  due  to  the  contact, 

r1  =  ^(S1-St1),  (5.7) 

which  represents  the  rate  at  which  carriers  are  injected  into  the  channel  for  each 
lattice  point  connected  to  the  leads.  The  broadening  matrices  become  important 
in  subsequent  sections  when  calculating  transport  dynamics  of  the  system. 


5.3  Retarded  Green’s  Function 


Given  these  definitions,  the  retarded  Green’s  function  can  be  computed  as 

G r(E)  =  [(E  +  iO+)I{NxNvNt)  -  H  -  SL  -  Sr]-1.  (5.8) 


This  quickly  becomes  intractable  for  large  system  sizes,  but  the  sparse,  block 
tridiagonal  structure  of  the  Hamiltonian  (see  Equations  5.2  and  5.3)  can  be  ex¬ 
ploited  to  convert  the  0((NxNyNz )3)  inversion  computation  to  0(NxNyNz)  using 
the  recursive  Green’s  function  (RGF)  method.  The  matrix  to  be  inverted  has  the 
following  block  tridiagonal  structure: 


G  r(E) 


Do  t0,i 

tl,0  Di  ti;2 

f  2,1  D2  t23 


f  AT— 2, AT— 3 


D 


N—2 


»N—l,N—2 


t>N-2,N-l 

Dat_i 


(5.9) 
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The  recursive  routine  involves  a  forward  and  backward  recursive  sweep.  For  the 
forward  sweep,  I  define  a  new  set  of  matrices, 


go.o  =  (Do)"1 

g if  =  (Dj  -  i  E  1,  2...N  -  1. 


(5.10) 


The  backward  recursion  yields  the  on-diagonal  blocks  of  the  retarded  Green’s 
function, 


GR  _  rrLK 

N—1,N—1  ~  SN-1,N-1 


LR 


Gg  =  +  t,,,+iG*  M+1t,+i,,g“), i  6  N  -  2, N  -  3...1, 0, 

where  I  is  the  identity  matrix  of  equivalent  size.  The  off-diagonal  blocks  are 


(5.11) 


Gf  J-«8W3ftw,  (5.12) 

’  ,Gf.  i>j 

The  entire  matrix  can  thus  be  resolved  with  the  recursive  algorithm.  The  blocks 
to  be  inverted  are  block  tridiagonal  themselves,  so  that  the  algorithm  can  be  used 
in  a  recursive  fashion  within  itself.  The  resulting  complexity  of  the  computation 

is  0(NxNyNz). 

The  electron  and  hole  correlations  are  given  by 


GniE)  =  GR[fLTL  +  fBTR}(GRy 
G P(E)  =  G*[(l  -  fL)VL  +  (1  -  fR)TR}( Gr)\ 


where  f^,  dependent  on  the  contact  chemical  potential,  ///,,  is  given  by  the  Fermi- 
Dirac  distribution, 


h(E) 


1 

e(E-iiL)/kT  _|_  ]_  ' 


(5.14) 
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5.4  Observables 


All  desirable  results  can  be  calculated  from  the  electron  and  hole  correlated  Green’s 
functions.  These  include  the  following: 

Tlr(E)  =  Tr(TLGBTR(GRy) 

2  roo 

JLR  =  jJ  TLR(E)[fL(E)  ~  fR(E)]dE  (5  15) 

/°°  df1 

T lr.  denotes  transmission  from  left  to  right  contacts,  and  current  I lR  is  the  integral 
over  the  transmission  for  which  the  Fermi-Dirac  distributions  of  the  contacts  are 
unequal.  The  spatially  resolved  electron  density,  phj,  is  an  integral  over  the  elec¬ 
tron  correlated  Green’s  function,  which  is  used  in  the  Poisson  solver  to  achieve 
a  self-consistent  solution.  Hole  density  is  the  same  integral  over  the  hole  corre¬ 
lated  Green’s  function,  so  that  the  sum  of  the  two  at  any  on-diagonal  point  is  the 
localized  density  of  states  for  the  particular  lattice/orbital. 

The  current  can  likewise  be  resolved  spatially,  so  that  the  total  current  from  any 
point  to  any  other  point  on  the  lattice  is 


ip  f  (i  f! 

Wi  =  j  /  -  Hj+hjG]J+1(E) 


(5.16) 


In  the  case  of  graphene,  the  Hamiltonian  is  only  nonzero  for  nearest  neighbors,  so 
j  +  1  denotes  a  lattice  point  that  is  only  one  point  away  from  j  in  any  direction. 
All  elements  in  Eq.  5.16  are  scalars  because  one  orbital  sufficiently  describes  each 
lattice  point.  The  total  current  between  lattice  points  is  thus  an  integral  of  scalar 
values  rather  than  an  integral  sum  over  matrices. 

As  a  verification,  the  sum  of  all  longitudinal  current  over  any  slice  in  the  y-z 
plane  is  constant  for  any  point  in  x  and  equivalent  to  Ilr.  All  matrix  elements  are 
verified  for  small  system  sizes  using  conventional  matrix  inversion/multiplication 
methods  to  ensure  the  validity  of  the  RGF  algorithm. 
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CHAPTER  6 


ANALYSIS  OF  TRANSPORT  IN 
SUPERFLUID  DOUBLE-LAYER 
GRAPHENE 


6.1  Approach 

The  system  of  interest  is  depicted  in  Fig.  4.3.  This  thesis  studies  two  electrostat¬ 
ically  doped  zigzag  graphene  monolayers  that  are  assumed  to  be  perfectly  regis¬ 
tered  with  one  another  and  separated  by  a  thin  dielectric.  Contacts  along  the  edges 
of  each  of  the  layers  inject  and  extract  current,  and  top  and  bottom  gates  manipu¬ 
late  the  quasiparticle  concentrations  in  each  of  the  layers.  In  this  way,  the  proper 
quasiparticle  concentrations  predicted  by  many -body  theory  [38]  may  be  tuned  to 
induce  a  superfluid  phase  transition.  The  top  and  bottom  gates  are  separated  from 
the  graphene  layers  by  20  nm  of  Si02.  The  two  graphene  layers  are  separated 
from  one  another  by  a  1  nm  S  i02  spacer  dielectric  to  be  in  the  predicted  regime  of 
superfluidity.  The  oxide  regions  are  assumed  to  be  perfect  in  the  sense  that  they  do 
not  contain  any  stray  charges  and  have  perfect  interfaces  with  the  graphene  layers. 
I  choose  the  x  direction  to  lie  along  the  length  of  the  system,  the  y  direction  to  lie 
along  the  width,  and  the  z  direction  along  the  depth.  Each  monolayer  is  30  nm 
long  by  10  nm  wide.  These  dimensions  were  selected  rather  than  larger  system 
sizes  so  that  the  atomistic  transport  properties  of  an  increased  number  of  disorder 
distributions  can  be  calculated  within  a  reasonable  amount  of  time. 

The  top  and  bottom  gates  (VTg  =  —Vbg)  arc  are  gated  to  effect  individual 
carrier  concentrations  of  1013  cm”2  in  each  layer,  which  corresponds  to  a  Fermi 
energy  in  the  top  (bottom)  layer  of  0.4  eV  (—0.4  eV).  The  gate  bias  conditions 
and  the  interlayer  separation  have  been  chosen  so  as  to  satisfy  the  conditions  for 
room  temperature  pseudospin  ferromagnetism  [38].  Here,  I  focus  on  the  regime 
where  the  layer  electron  and  hole  populations  place  the  system  firmly  in  the  dense 
electron-hole  regime  where  the  electron-hole  pairs  are  expected  to  form  a  BCS- 
type  state.  In  practice,  this  is  a  nontrivial  density  to  achieve,  but  is  necessary 
for  condensation.  The  electric  field  created  (8  x  106  V/cm)  between  the  top  and 
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bottom  monolayers  is  slightly  below  the  dielectric  breakdown  voltage  for  Si02 
and  is  theoretically  feasible.  The  transport  characteristics  of  such  a  system  are 
in  stark  contrast  to  the  dilute  limit  of  electron-hole  densities,  which  can  be  de¬ 
scribed  as  a  weakly  interacting  Bose  system  of  excitons  [41].  In  order  to  collect  a 
sufficient  statistical  distribution  on  the  effects  of  vacancy  distribution  on  the  inter¬ 
layer  transport  properties  of  double-layer  graphene,  several  randomized  vacancy 
configurations  for  each  concentration  are  examined  in  this  work. 

The  three-dimensional  Hamiltonian,  as  described  in  Chapter  5,  begins  with  the 
atomistic  tight-binding  description  of  an  individual  graphene  monolayer, 

Htl  =  ^  T  I  *)0’  I  +v*  I  *>(*  l>  (6-1) 

<i,j> 

where  lattice  points  i  and  j  are  first  nearest  neighbors,  r  =  —3.03  eV  is  the  near¬ 
est  neighbor  hopping  energy  for  the  pz  orbital  of  graphene,  which  allows  for  the 
unique  low-energy  linear  dispersion  at  the  K  and  K'  points  in  the  Brillouin  zone 
shown  in  Fig.  4.2.  Hopping  among  further  nearest  neighbors  and  other  orbitals  is 
neglected,  because  nearest  neighbor  p2  orbital  hopping  is  the  predominant  inter¬ 
action  for  graphene  in  the  probed  energy  range  of  E#  ±  1  eV.  The  on-site  potential 
energy  Vt  =  0(f))  is  calculated  via  a  three-dimensional  Poisson  solver  given  in 
Eq.  5.1.  The  contact  self-energy  matrices  use  a  phenomenological  model  to  simu¬ 
late  a  generic  metal  contact  with  a  constant  density  of  states  [30;  42].  This  model 
captures  the  basic  self-energy  needed  to  appropriately  simulate  a  metal  contact 
without  taking  into  account  multiple  orbitals  or  complex  interface  problems  such 
as  lattice  mismatch  or  Schottky  barrier  height. 

The  monolayer  Hamiltonian  is  expaned  to  the  double-layer  Hamiltonian  by 
coupling  the  top  and  bottom  monolayers  with  the  following  Bogoliubov-de  Gennes 
(BdG)  Hamiltonian, 

+  ^  jj,  •  A  <g>  crM,  (6.2) 

l_i=x,y,z 

with  the  interlayer  interactions  including  both  single  particle  tunneling  and  the 
mean-field  many-body  contribution,  A,  coupling  the  two  layers  using  a  local  den¬ 
sity  approximation.  In  Eq.  6.2,  //  represents  a  vector  that  isolates  each  of  the 
Cartesian  components  of  the  pairing  vector,  represents  the  Pauli  spin  matrices 
in  each  of  the  three  spatial  directions,  and  g)  represents  the  Kronecker  product. 


T-t-Bdc  — 


Htl  0 
0  Hbl 
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In  order  to  correctly  account  for  the  dynamics  of  the  double-layer  graphene 
system,  I  first  include  the  many-body  interlayer  interactions  in  the  Hamiltonian. 
Within  the  Hartree-Fock  mean-field  approximation,  the  interlayer  interactions  may 
be  defined  through  the  expectation  value  of  the  full  Hamiltonian, 


(til  Hb<ig  |tj)  UmexcSij. 


(6.3) 


The  graphene  monolayers  are  assumed  to  be  perfectly  registered,  so  that  electrons 
at  site  i  in  the  top  layer  ((til)  only  bind  with  holes  at  site  j  on  the  bottom  layer 
(ltj))  when  i  =  j.  U  is  the  strength  of  the  interlayer  on-site  Coulomb  interaction, 
whose  selected  value  is  addressed  later  in  this  section. 

The  magnitude  of  the  order  parameter  resulting  from  the  analysis  of  BECs  of 
indirectly  bound  excitons  is  rnexc.  It  is  proportional  to  the  off-diagonal  terms  in 
the  single-particle  density  matrix,  which  is  represented  as  [13;  39] 


P  = 


Pn  Pn 
Pi t  Ait 


(6.4) 


The  on-diagonal  density  matrix  (ptt>Pw)  corresponds  to  the  associated  electron 
and  hole  densities  of  the  top  and  bottom  monolayers,  respectively.  The  order 
parameter  mexc  can  now  be  defined  as  a  function  of  the  interlayer  component  of 
the  density  matrix, 


mexc  =  Pn  +  Pi  t  =  2Re(Pn ), 

™vexc  =  -*Au  +  iPi  t  =  2  Im(pn). 


(6.5) 


The  density  matrix  is  directly  calculated  within  the  NEGF  formalism.  It  is  con¬ 
sequently  both  an  input  and  an  output  of  the  simulation,  and  the  algorithm  must 
iterate  over  the  above  mean-field  equations,  in  conjunction  with  the  Poisson  equa¬ 
tion  for  electrostatics,  to  obtain  a  self-consistent  solution  with  compatible  particle 
densities  and  potential  profile.  The  Broyden  method  [43]  is  utilized  to  expedite 
numerical  convergence. 

Equation  6.2  can  be  expanded  to  show  a  simplified  form  for  the  BdG  Hamil¬ 
tonian  in  which  interlayer  interactions  are  expressed  in  terms  of  their  directional 
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components  [44], 


T^BdG  ~ 


H- 


A.  At  -  *A, 


TL  i  x  <"-» y 

Ax  +  iAy  Hbl  —  A z 


(6.6) 


The  directional  components  of  the  interlayer  interactions  A  are  expressed  as 


A*  =  (Asas  +  Umxexc) 

Ay  =  Umyexc  (6.7) 

|  i)(i  |  -Vj.  |  i)(i  |). 

The  on-diagonal  term  in  the  interlayer  interactions,  Az,  is  due  to  screening  caused 
by  the  unbound  carriers  in  each  monolayer,  and  acts  to  separate  the  two  Fermi 
surfaces.  The  value  of  the  single  particle  tunneling  energy,  Asas,  is  proportional 
to  the  probability  of  a  single  electron  tunneling  through  the  thin  dielectric  and 
recombining  with  a  hole.  Single-particle  tunneling  is  an  adverse  event;  thus  it  is 
desirable  to  have  a  very  thick  barrier  in  between  the  two  graphene  layers  to  max¬ 
imize  the  lifetime  of  the  indirectly  bound  excitons.  However,  a  strong  interlayer 
Coulomb  interaction  is  needed  to  drive  the  superfluid  phase  transition,  which  ne¬ 
cessitates  a  thin  barrier.  In  this  work,  I  set  Asas  =  1  //cV  sufficiently  small  so  that 
the  lifetime  of  the  indirectly  bound  exciton  is  long  enough  to  observe  condensa¬ 
tion  but  not  so  small  as  to  require  an  intractably  large  number  of  iterations  before 
self-consistency  is  reached. 

Results  of  previous  many-body  calculations  show  the  value  of  the  order  pa¬ 
rameter  is  approximately  one-tenth  the  Fermi  energy  [38].  The  simulations  show 
that  this  value  of  mexc  corresponds  to  an  interlayer  coupling  strength  of  2.0  eV. 

Although  this  is  less  than  the  unscreened  mean-field  interaction  strength  [13] 

2 

(U  =  2-;  ~  4.6  eV),  it  is  more  plausible  because  it  factors  in  damping  effects 
due  to  screening  [45].  The  interlayer  Coulomb  interaction  is  thus  set  to  U  =  2.0 
eV  in  all  simulations  to  be  compatible  with  the  mean-field  theory.  The  magnitude 
of  the  order  parameter,  using  the  expectation  values  of  the  density  matrix  in  Eq. 
6.5,  is 

\mexc\  =  \/  {mexc)2  +  (m exc )2  (6.8) 

The  phase  of  the  order  parameter  is  similarly  dependent  on  the  same  expectation 
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values  of  the  density  matrix  [44], 


Yexc 


tan  1 


m 

m 


y 

exc 

x 

exc _ 


(6.9) 


As  the  interlayer  transport  properties  are  similar  to  the  case  of  a  Josephson  junc¬ 
tion,  the  phase  is  expected  to  be  <f>exc  =  7r/2  at  zero  bias  to  maximize  interlayer 
current  [46;  47].  The  quasiparticle  and  condensate  current  densities,  the  artifacts 
of  the  Andreev  reflection  action,  are  proportional  to  the  spatial  phase  gradient  and 
magnitude  of  the  order  parameter  [13]. 


6.2  The  Ideal  System 

All  simulations  are  run  in  the  drag-counterflow  geometry  ( Vtl  =  — Vtr ;  Vbl  = 
Vbr  =  0  V)  [13]  to  induce  an  effective  current  through  the  bottom  layer  only 
when  the  condensate  exists,  as  described  in  Section  3.1.  All  non-equilibrium  con¬ 
figurations  described  below  are  in  this  geometry,  with  the  left  and  right  contacts 
on  the  top  monolayer  set  to  the  magnitude  of  the  counterflow  bias  parameter.  The 
entire  process  results  in  a  condensate  current,  due  to  the  propagating  exciton,  and 
a  quasiparticle  current,  caused  by  the  injected  and  retro-reflected  individual  car¬ 
riers.  The  quasiparticle  current  is  only  nonzero  within  the  coherence  length  ( Lc ), 
the  maximum  length  an  injected  particle  penetrates  into  the  superfluid  gap  before 
triggering  the  exciton  scattering  event.  The  microscopic  attributes  of  the  super¬ 
fluid  and  the  dynamics  of  the  Andreev  reflection  may  not  be  readily  apparent  in 
experiments,  but  will  have  a  macroscopic  effect  on  the  observable  interlayer  cur¬ 
rent. 

These  observable  interlayer  and  intralayer  transport  properties  are  presented  for 
an  ideal  system  in  Fig.  6.1.  The  currents  are  odd  functions  of  the  bias,  reflecting 
the  system’s  ambipolar  nature.  At  low  bias,  interlayer  current  is  linear  and  inter¬ 
layer  conductance  is  constant,  because  injected  carriers  see  the  superfluid  gap  and 
trigger  the  Andreev  event  discussed  above.  Low-energy  injected  carriers  are  un¬ 
able  to  pass  through  the  superfluid  gap,  and  intralayer  current  is  negligible  within 
this  range.  States  above  the  superfluid  band  gap  do  begin  to  form  at  \E\  >  0.14 
eV,  as  seen  in  the  superfluid  band  structure  in  Fig.  6.2a.  Interlayer  transmission 
quickly  vanishes  when  higher-energy  states  become  accessible  in  the  spectrum; 
transport  along  the  monolayer  dominates  in  these  ranges.  This  is  evidenced  by  the 
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Figure  6.1:  Interlayer  and  intralayer  current  plots  for  the  ideal  system  as  a  function 
of  counterflow  bias  (solid  black  line  and  red  dashed  line,  respectively). 


emergence  of  intralayer  current  and  a  drop  in  interlayer  conductance  beginning  at 
a  counterflow  bias  of  0.14  V. 

Beyond  the  critical  counterflow  bias,  0.22  ±  0.01  V  in  the  ideal  system,  the  su¬ 
perfluid  can  no  longer  adjust  its  phase  to  accommodate  the  current  flow.  When 
self-consistency  is  lost,  only  single -particle  tunneling  contributes  to  interlayer 
tunneling.  The  small  magnitude  of  the  resultant  interlayer  current,  less  than  1 
pA  for  a  bias  of  VTl  =  —  VTr  =  0.25  V,  is  negligible  compared  to  the  many- 
body  contribution  to  the  interlayer  current  observed  when  the  system  obtains  self- 
consistency.  Additionally,  past  the  critical  current,  transport  characteristics  be¬ 
come  time-dependent  [48],  and  the  mean  field  equations  no  longer  represent  a 
valid  description  of  the  physics.  Although  a  transient  remnant  of  the  superfluid 
may  yet  persist  beyond  the  critical  bias,  the  many-body  portion  of  the  Hamilto¬ 
nian  is  approximately  zero  and  the  time  dependent  oscillation  in  the  interlayer 
transport  relationship  cannot  be  captured  in  the  steady-state  simulations. 

The  phase  transition  is  observed  in  Fig.  6.2b,  which  shows  the  dispersion  for 
the  ideal  system  after  critical  current  has  been  passed.  There  is  no  appreciable 
change  to  the  band  structure  within  each  phase.  When  no  superfluid  exists,  each 
monolayer  has  a  vanishing  band  gap  and  linear  dispersion,  as  seen  in  Fig.  4.2.  As 
expected,  the  Dirac  points  occur  in  the  normal  phase  at  the  K’  and  K  points  in  the 
Brillouin  zone  with  a  Fermi  energy  of  ±0.4  eV.  The  band  gap  vanishes  at  biases 
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k  ■  a/71 


Figure  6.2:  The  dispersion  relation  for  an  ideal  double-layer  graphene  system  in 
the  transport  direction  (a)  at  a  counterflow  bias  of  VTL  =  —VTr  =  0.05  V,  which 
is  below  the  critical  current  and  (b)  at  a  counterflow  bias  of  Vtl  =  —Vtr  =  0.23 
V,  which  is  above  the  critical  current. 
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larger  than  this  critical  transition  value.  The  closing  of  the  band  gap  after  the 
phase  transition  allows  for  more  low-energy  states,  causing  a  spike  in  intralayer 
current  well  beyond  the  range  of  interlayer  transport  in  the  superfluid  phase. 

The  band  structure  in  the  transport  direction  of  the  system  is  derived  [49]  from 
a  small  portion  of  the  converged  Hamiltonian  by  assuming  mexc  is  periodic  in 
the  transport  direction.  This  holds  true  only  for  the  ideal  system,  where  mexc  is 
smooth  and  consistent  throughout  the  channel  with  a  magnitude  near  10%  of  the 
Fermi  energy  and  a  phase  of  tt/2.  Randomly  placed  impurities  and  vacancies, 
however,  break  translational  symmetry  so  that  the  Hamiltonian  is  no  longer  pe¬ 
riodic  in  the  transport  direction.  The  dispersion  relations  of  disordered  systems 
thus  cannot  be  calculated. 

The  ideal  system’s  critical  current  is  calculated  to  be  Ic  =  ±18.4  //A,  occurring 
at  a  counterflow  bias  of  ±0.22  V.  The  analytic  approximation  for  the  critical  cur¬ 
rent  in  a  system  where  coherence  length  is  smaller  than  system  length  [48]  can  be 
evaluated  by  ensuring  the  condensate  satisfies  the  elliptic  sine-Gordon  equation, 

A2V2</>  -  sin{<j>)  =  0.  (6.10) 


When  this  equation  is  solved  in  the  static  (DC)  case,  the  critical  current  of  the 
system  simplifies  to  a  relatively  simple  expression, 


h 


r^j 


eWmexc 

hLc 


(6.11) 


Therefore,  for  a  system  of  width  W  =  10  nm  with  a  coherence  length  of  Lc  «  5 
nm  and  order  parameter  magnitude  of  ps  ~  0.04  eV,  the  analytic  critical  current 
is  roughly  1^  «  19.5  pA.  This  is  in  good  agreement  with  the  simulated  critical 
current  for  an  ideal  system.  The  reduction  in  the  critical  current  is  expected  in 
the  system  as  reflections  off  the  superfluid  gap  lead  to  enhanced  differences  in  lo¬ 
cal  interlayer  interaction  terms  which  enter  into  the  Hamiltonian  self-consistently 
through  the  A,  term  and  serve  to  energetically  separate  the  layers  [44]. 

Figure  6.3  plots  the  interlayer  quasiparticle  current  densities  along  the  transport 
direction  averaged  along  the  width  of  the  system  in  an  ideal  system  with  a  coun¬ 
terflow  bias  of  0.05  V  and  0.20  V,  just  before  the  phase  transition  to  a  normal  state. 
Quasiparticle  current  tunneling,  which  is  part  of  the  same  process  that  launches 
the  exciton  through  the  channel,  is  highest  near  the  contacts  where  carriers  are 
injected,  and  evanescently  decays  into  the  channel,  so  that  carriers  which  provoke 
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Length  (nm) 

Figure  6.3:  The  interlayer  quasiparticle  current  density  along  the  transport  direc¬ 
tion  in  the  ideal  system  averaged  along  the  width  of  the  system  for  counterflow 
biases  of  0.05  V  (red  dashed  line)  and  0.20  V  (blue  solid  line).  The  current  density 
is  only  nonzero  within  the  coherence  length,  which  by  inspection  is  approximately 
5  nm.  Quasiparticle  current  density  is  approximately  odd  about  the  center  of  the 
system  at  low  bias,  but  symmetry  is  lost  when  a  high  contact  bias  limits  carrier 
concentration  on  one  side  of  the  system. 
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interlayer  transmission  are  only  present  a  small  distance  into  the  channel.  This 
distance,  termed  the  coherence  length  ( Lc ),  is  roughly  5  nm  by  inspection  of  the 
plot. 

Note  in  Fig.  6.3  that  interlayer  quasiparticle  current  magnitude  is  nearly  equiv¬ 
alent  at  low  bias  but  becomes  very  asymmetric  near  the  critical  transition  point. 
The  Poisson  equation  in  the  counterflow  bias  configuration  restricts  the  magni¬ 
tude  of  the  order  parameter  because  the  negatively  biased  top  right  contact  locally 
decreases  the  density  of  electrons  able  to  pair  into  the  condensate.  Conversely, 
applying  a  bias  across  the  hole-doped  bottom  layer  will  limit  superfluid  density 
and  quasiparticle  current  near  the  left  (positively  biased)  contact.  The  decreased 
density  results  in  a  smaller  interlayer  current  density  on  one  side,  as  seen  in  the 
case  with  higher  counterflow  bias.  All  interlayer  currents  are  the  currents  gener¬ 
ated  by  the  contacts  on  the  right  side  of  the  system,  because  this  is  the  side  with 
limited  exciton  density  is  the  transport  bottleneck.  When  the  disparity  in  conden¬ 
sate  currents  from  each  contact  reaches  a  critical  limit,  steady-state  supercurrents 
are  no  longer  possible  [13]. 


6.3  Disordered  Systems 

I  now  seek  to  understand  how  the  behavior  of  the  interlayer  transport  properties 
change  as  disorder  is  introduced  into  the  system.  The  main  result  of  the  analysis 
is  shown  in  Fig.  6.4,  which  plots  the  statistical  variation  in  critical  current  as  a 
function  of  vacancy  concentration.  Two  distinct  categories  of  disorder  in  the  top 
layer  are  considered.  The  critical  current  obtained  for  the  ideal  case  is  plotted  in 
the  middle  of  the  figure,  at  a  vacancy  concentration  of  0%.  To  the  left  of  the  ideal 
case,  the  statistical  variation  in  the  critical  current  is  plotted  for  random  vacancy 
concentrations  which  are  beyond  the  coherence  length,  Lc,  from  either  contact 
on  the  top  layer.  In  this  case,  an  increase  in  the  random  vacancy  concentrations 
causes  a  decrease  in  the  critical  current  approximately  following  a  square  root  de¬ 
pendence.  Furthermore,  the  mean  of  the  statistical  distribution  of  critical  currents 
associated  with  various  top  layer  vacancy  distributions  varies  only  by  approxi¬ 
mately  5%  even  as  the  concentration  of  vacancies  is  increased. 

The  statistical  change  of  the  critical  current  associated  with  disorder  within  a 
distance  of  Lc  from  either  top  layer  contact  is  plotted  on  the  right  hand  side  of 
Fig.  6.4.  Here  quite  different  behavior  is  observed  compared  to  the  situation 
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Figure  6.4:  Statistical  variation  of  critical  interlayer  current  for  several  runs  at 
the  bias  just  before  the  condensate  is  broken.  Vacancies  which  are  more  than  a 
distance  of  one  coherence  length  away  from  the  top  layer  contacts  are  to  the  left  of 
the  ideal  (0%)  case,  depicted  by  a  black  circle,  while  vacancies  which  are  within 
one  coherence  length  from  the  top  layer  contacts  are  to  the  right.  The  excitonic 
condensate  is  lost  in  both  cases  for  top  layer  vacancy  concentrations  larger  than 
4%.  The  solid  lines  represent  analytic  calculations  using  Eq.  6.11  with  modified 
values  obtained  from  calculations  as  described  in  Section  6.3. 
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Figure  6.5:  Set  of  I-V  curves  for  systems  with  varying  concentrations  of  randomly 
placed  vacancies  at  least  5  nm  away  from  the  top  layer  contacts. 

where  the  vacancies  are  further  than  Lc  from  the  top  layer  contacts.  Now  the 
critical  current  has  a  very  distinct  linear  decrease  with  respect  to  an  increase  in 
vacancy  concentrations.  There  are  very  large  variations  in  the  location  of  the 
critical  current,  which  is  dependent  on  the  locations  of  the  top  layer  vacancies.  In 
the  subsequent  sections,  I  will  explore  the  physics  of  these  two  distinct  regions 
and  explain  the  observed  dependence  of  the  critical  current  seen  in  each  situation. 


6.3.1  Layer  Disorder  beyond  Lc 

Disorder  is  modeled  by  randomly  placing  vacancies  within  the  top  layer  of  the 
double-layer  system.  Vacancies  are  modeled  by  setting  all  hopping  to  a  missing 
atom  to  zero  in  the  model  Hamiltonian.  This  effectively  blocks  any  interaction 
with  the  vacancy  by  setting  the  tight-binding  overlap  integral  of  the  spatial  pz  or¬ 
bital  states  to  zero  [50].  A  fixed  percentage  of  carbon  atoms  are  randomly  removed 
from  the  specified  region  of  the  top  monolayer,  leaving  the  bottom  monolayer  un¬ 
perturbed,  and  the  transport  calculations  from  the  ideal  case  are  repeated. 

The  interlayer  current  is  plotted  as  a  function  of  voltage  bias  in  Fig.  6.5  for  var¬ 
ious  concentrations  of  vacancies  in  the  top  layer  along  with  the  ideal  double-layer 
case  for  comparison.  At  low  bias,  there  is  little  reduction  in  the  interlayer  conduc¬ 
tivity  regardless  of  vacancy  concentration  or  vacancy  location  past  the  coherence 
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Figure  6.6:  Magnitude  of  the  order  parameter  for  a  system  with  1%  vacancies  in 
the  channel.  The  condensate  is  very  close  to  that  of  the  ideal  case,  with  a  drop 
in  magnitude  only  occuring  locally  near  disorder.  As  expected,  \mexc\  «  0.1  EF 
away  from  disorder. 


length.  The  phase  transition  to  two  normal,  incoherent  Fermi  liquids,  however, 
occurs  at  a  lower  bias  than  in  the  ideal  case.  Channel  disorder  decreases  critical 
current  by  20%  for  1%  vacancy  concentration  to  30%  for  vacancy  concentrations 
of  4%.  Beyond  top  layer  vacancy  concentrations  of  4%,  no  self-consistent  solu¬ 
tion  is  found  and  the  interlayer  transport  is  dominated  by  single-particle  tunneling 
events.  Resultant  interlayer  currents  are  orders  of  magnitude  smaller  because  the 
majority  of  the  current  injected  into  the  system  now  flows  across  the  graphene 
layers. 

In  order  to  better  understand  these  reductions  in  interlayer  current  as  the  top 
layer  vacancy  concentration  is  increased,  I  examine  the  magnitude  of  the  order 
parameter,  \mexc\.  In  Fig.  6.6,  the  magnitude  of  the  order  parameter  is  plotted 
for  a  random  vacancy  concentration  of  1%.  In  this  situation,  the  order  parameter 
magnitude  does  not  remain  constant  over  the  entire  system  at  10%  EF,  as  was  the 
case  in  the  ideal  system.  The  vacancies  locally  destroy  the  condensate  [45]  and 
reduce  \mexc\  at  surrounding  points  up  to  a  distance  of  0.5  nm,  or  three  to  four 
nearest  neighbors  from  the  vacancies.  However,  no  long-range  effect  is  seen  when 
vacancies  are  isolated  from  one  another  by  more  than  approximately  2  nm.  At  this 
particular  vacancy  concentration,  mexc  is  reduced  by  40%  over  an  appreciable  area 
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Figure  6.7:  Magnitude  of  the  order  parameter  for  a  system  with  4%  vacancies  in 
the  channel.  Magnitude  clearly  drops  nonlinearly  with  the  increase  in  disorder, 
but  quasiparticle  current  is  still  able  to  form  near  the  contacts. 


of  the  system  as  a  result  of  these  vacancies. 

This  situation  is  to  be  contrasted  with  Fig.  6.7,  in  which  mexc  is  plotted  for 
a  random  top  layer  vacancy  concentration  of  4%.  In  this  case,  significant  areas 
clearly  emerge  where  concentrated  disorder  has  long-range  effects  on  superfluid 
density.  These  areas  of  high  vacancy  concentrations  in  Fig.  6.7  give  rise  to  values 
for  \mexc\  that  are  less  than  20%  the  ideal  value  over  significant  areas  of  the  sys¬ 
tem.  Superfluid  magnitude  remains  at  its  ideal  level  near  the  contacts,  sufficiently 
quarantined  from  the  vacancies. 

Nevertheless,  the  root  cause  of  the  sublinear  behavior  in  the  critical  current  with 
respect  to  an  increase  in  the  top  layer  vacancy  concentration  remains  unresolved, 
from  Eq.  6.11.  The  localized  density  of  states  (LDOS)  [51]  is  examined  to  ex¬ 
plain  the  sublinear  behavior.  Vacancies  induce  a  LDOS  similar  to  the  case  of  an 
impurity  [52;  53]  in  graphene  layers.  Figures  6.8  and  6.9  plot  the  LDOS  for  the 
ideal  case  and  for  the  case  of  4%  top  layer  vacancy  concentration  as  an  average 
along  the  disordered  area  of  each  monolayer.  The  average  shows  the  magnitude  of 
uncoupled  states  that  form  not  only  at  vacancies,  but  at  neighboring  locations  and 
in  the  ideal  bottom  layer  as  a  result  of  electrostatics.  A  stark  contrast  is  seen  in 
the  low-energy  LDOS  in  the  superfluid  phase  when  vacancies  are  introduced.  The 
LDOS  closely  resembles  those  of  monolayer  graphene  at  higher  energies  [54], 
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Figure  6.8:  Local  density  of  states  at  low  energy  for  an  ideal  system.  The  solid 
(dashed)  line  represents  the  top  (bottom)  layer.  No  localized  states  exist  below  a 
magnitude  of  0.14  eV  for  the  superfluid. 
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Figure  6.9:  Local  density  of  states  at  low  energy  for  a  system  with  4%  vacancies 
in  the  top  layer.  The  solid  (dashed)  line  represents  the  top  (bottom)  layer.  Low- 
energy  states  due  to  the  vacancies  introduce  a  weak  but  apparent  LDOS  in  the 
ideal  bottom  layer  as  well,  due  to  the  thinness  of  the  spacer  dielectric. 
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with  peaks  at  E  =  EF  ±  r.  The  ideal  top  and  bottom  monolayers  exhibit  per¬ 
fectly  antisymmetric  densities  of  states  so  that  equivalent  carrier  concentrations 
arise  in  the  oppositely  gated  top  and  bottom  monolayers.  No  states  exist  in  the 
superfluid  gap  in  the  ideal  double-layer,  as  is  expected  for  a  condensate  in  which 
all  of  the  quasiparticles  participate. 

Disorder  changes  the  transport  properties  of  the  system  by  introducing  mid¬ 
gap  states  in  the  top  layer,  as  seen  in  Fig.  6.9.  Localized  states  also  arise,  to 
a  lesser  degree,  in  the  bottom  layer  due  to  charge  pileup  induced  by  Coulomb 
attraction  through  the  thin  spacer  dielectric.  The  biased  top  and  bottom  gates, 
necessary  to  generate  the  sufficient  carrier  concentrations,  cause  the  undesirable 
occupation  of  the  mid-gap  states  up  to  the  Fermi  energy.  As  the  electrons  and 
holes  make  their  way  across  the  layers,  they  scatter  off  localized  states,  changing 
the  interlayer  phase  relationship  between  the  top  and  bottom  layers.  This  effect 
manifests  itself  as  a  change  in  the  Ay  term  in  Eq.  6.5.  As  a  result,  when  the 
decreased  interlayer  phase  relationship  is  input  into  the  calculation  of  mexc  in 
Eq.  6.8,  mexc  drops  from  its  ideal  value  of  0.041  eV  to  0.029  ±  0.002  eV  at  4% 
top  layer  vacancy  concentration,  a  decrease  of  roughly  30%.  The  result  is  thus 
that  the  reduced  interlayer  phase  component  decreases  mexc  in  a  rough  square 
root  correlation.  The  square  root  dependence  is  not  exact  and  saturates  at  higher 
disorder  concentrations  because  mexc  is  not  affected  near  the  contact  regions  and 
is  allowed  to  remain  at  its  ideal  value,  as  seen  in  Fig.  6.7.  When  the  average  value 
of  the  change  in  the  coupling  strength  is  added  into  the  calculation  of  the  critical 
current  in  Eq.  6.11  and  plotted  along  with  the  numerical  results  on  the  left  hand 
side  of  Fig.  6.4,  good  agreement  is  found  between  the  two  values.  The  average 
coupling  strength  does  not  change  appreciably  between  disorder  configurations, 
resulting  in  critical  current  approximations  with  small  deviations  relative  to  the 
contact  disordered  scenario. 


6.3.2  Layer  Disorder  within  Lc 

While  the  double-layer  graphene  system  is  moderately  robust  to  vacancies  deep  in 
the  channel,  this  is  not  the  case  when  the  top  layer  vacancies  occur  within  one  Lc 
of  the  contacts  on  the  top  layer.  The  interlayer  current  is  plotted  as  a  function  of  the 
random  top  layer  vacancy  concentration  in  Fig.  6.10,  where  the  vacancies  occur 
within  one  Lc  of  the  contacts  on  the  top  layer.  Vacancies  near  the  contacts  cause  a 
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Figure  6.10:  I-V  curve  for  systems  with  varying  concentrations  of  randomly 
placed  vacancies  within  5nm  of  the  top  layer  contacts. 


significant  drop  in  both  conductivity  and  the  voltage  bias  at  which  self-consistency 
is  lost.  Quasiparticle  current  density  is  now  dependent  on  the  location  of  disorder 
within  the  coherence  length,  because  magnitude  evanescently  decays  from  the 
contact  interface.  Interlayer  current  is  thus  sensitive  to  the  specific  location  and 
amount  of  disorder  within  the  coherence  length. 

Figure  6.11  shows  the  magnitude  of  the  order  parameter,  mexc,  for  a  case  of 
1%  vacancies  randomly  distributed  near  the  contacts.  Quasiparticle  current  den¬ 
sity  drops  with  the  presence  of  contact  disorder,  because  it  is  proportional  to  the 
magnitude  of  the  order  parameter.  Clearly,  groups  of  vacancies  appear  on  the  left 
and  right  side  of  the  contact  that  begin  to  show  a  nonlinear,  long-range  effect  on 
superfluid  density.  Top  layer  vacancies  closest  to  the  contact,  where  quasiparticle 
tunneling  magnitude  is  largest,  cause  the  biggest  detriment  to  the  magnitude  of  in¬ 
terlayer  current.  Despite  the  vacancies  reducing  the  space  in  which  quasiparticles 
may  be  injected  without  scattering,  the  condensate  is  able  to  form;  sites  within 
three  nearest  neighbors  of  a  single  defect  are  not  appreciably  affected.  Little  ran¬ 
domized  bunching  occurs  in  the  1%  case,  because  the  disorder  is  too  sparse  to 
create  significantly  different  scenarios.  Interlayer  current  remains  rather  robust, 
roughly  30%  smaller  than  ideal.  The  random  nature  of  the  placement,  however, 
causes  a  high  variance  in  interlayer  current  in  this  scenario,  which  is  discussed 
below. 
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Figure  6.11:  Magnitude  of  the  order  parameter,  in  eV,  for  a  system  subjected  to 
1%  vacancies  near  the  contacts.  A  single  defect  does  not  affect  the  condensate 
appreciably,  but  several  defects  nearby  can  drop  magnitude  significantly. 


Larger  deviations  from  the  ideal  interlayer  critical  current  are  found  as  the  top 
layer  vacancy  concentration  is  increased.  A  large  amount  of  clustering  occurs, 
which  leads  to  a  significant  reduction  in  the  available  space  into  which  quasiparti¬ 
cles  may  be  injected,  as  seen  in  Fig.  6.12,  and  generates  a  very  different  transport 
relationship  relative  to  channel  disorder.  Carriers  injected  into  the  system  see  a 
significant  reduction  in  the  area  in  which  quasiparticle  tunneling  can  occur  near 
the  contact,  significantly  reducing  conductivity.  There  is  a  linear  dependence  on 
the  width  of  the  system  in  Eq.  6.1 1,  and  this  provides  a  very  simple  explanation 
for  the  physics  of  the  linear  decrease  in  the  critical  current  observed  in  Fig.  6.4. 
This  information  allows  us  to  conclude  that  graphene  with  vacancies  within  the 
coherence  length  effectively  reduces  the  width  over  which  quasiparticle  tunneling 
occurs,  and  causes  a  linear  decay  in  Ic  with  respect  to  contact  disorder  strength 
until  the  condensate  vanishes. 

Figure  6.13  shows  how  the  average  interlayer  quasiparticle  current  density  qual¬ 
itatively  shifts  when  top  layer  vacancies  are  included  within  the  coherence  length 
of  the  top  layer  contacts.  Equivalent  disorder  concentrations  generate  disparate 
condensate  currents  on  each  side  because  vacancies  are  randomly  configured  to  be 
more  closely  lumped  near  one  contact  than  the  other.  Whereas  superfluid  excitons 
are  able  to  permeate  disorder  as  long  as  the  condensate  exists,  an  increasing  pro- 
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Figure  6.12:  Magnitude  of  the  order  parameter,  in  eV,  for  a  system  subjected  to 
4%  vacancies  near  the  contacts.  The  contacts  are  too  saturated  with  defects  for  a 
large  amount  of  tunneling  to  occur,  and  interlayer  conductivity  suffers. 
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Figure  6.13:  Quasiparticle  current  density  for  a  system  with  4%  contact  disorder 
with  a  counterflow  bias  of  0.04  V,  just  before  the  phase  transition.  No  condensate 
exists  near  the  contacts,  so  the  largely  degraded  quasiparticle  tunneling  that  does 
occur  only  happens  beyond  the  disorder. 


40 


portion  of  bare  electrons  and  holes  cannot  penetrate  the  contact  disorder  to  reach 
the  condensate.  The  critical  bias  at  which  the  phase  transition  occurs  decreases 
as  disorder  increases  because  of  the  discrepancy  in  condensate  current  from  each 
contact.  The  on-site  potential  is  greater  at  vacancies  near  the  contacts  due  to  VTR 
and  Vtl •  The  LDOS  is  thus  more  quickly  occupied,  and  further  accelerates  the 
phase  transition.  As  a  result,  critical  current  drops  approximately  by  an  average 
of  30%,  45%,  60%,  and  80%  for  vacancy  concentrations  of  1%,  2%,  3%,  and  4%, 
respectively.  The  accuracy  of  the  conclusion  can  be  seen  on  the  right  hand  side  of 
Fig.  6.4  where  the  data  is  superimposed  with  a  model  of  the  reduction  of  the  crit¬ 
ical  current  expected  when  each  vacancy  is  assumed  to  remove  0.5  nm  from  the 
effective  width  of  the  system  in  the  affected  region  for  a  given  vacancy  concen¬ 
tration.  The  theory  and  numerical  calculations  prove  to  be  in  good  agreement;  the 
model  falls  within  the  error  bars  of  the  simulations,  where  deviation  is  the  result 
of  non-linear  randomized  bunching  of  disorder. 

It  is  clear  in  Fig.  6.4  that  a  large  variance  in  critical  current  exists  when  va¬ 
cancies  lie  within  the  coherence  length  of  the  contacts  in  the  top  layer.  This  is 
due  to  the  fact  that  the  quasiparticle  current  flow,  in  this  system,  is  concentrated 
close  to  the  contacts  in  the  middle  of  the  channel,  and  evanescently  decays  into 
the  system.  This  is  best  seen  in  Fig.  6.11  where  the  magnitude  of  the  order  param¬ 
eter  is  significantly  reduced  from  the  bulk  value  near  the  edge  of  the  system  both 
near  and  away  from  the  contacts.  Moreover,  vacancies  that  are  closely  lumped 
together  exacerbate  the  scattering  to  an  even  greater  degree  because  vacancies  in 
close  proximity  to  one  another  have  a  significant  non-local  effect  on  the  perturbed 
mexc,  as  seen  in  Fig.  6.12.  Also  in  Fig.  6.12,  the  variance  in  critical  current  decays 
at  4%  vacancies  because  quasiparticle  tunneling  is  almost  completely  prevented 
at  such  a  high  concentration  of  vacancies.  As  in  the  channel-disordered  case,  no 
self-consistent  solution  was  found  for  vacancy  concentrations  greater  than  4%. 
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CHAPTER  7 


CONCLUSION 


Superfluid  transport  is  an  observed  phenomenon  in  BECs  in  which  bosons  flow 
through  a  channel  without  dissipation.  Many-body  transport  in  such  a  system 
can  lead  to  efficiencies  and  speeds  beyond  the  fundamental  limit  of  conventional 
CMOS  technology.  The  performance  enhancement  could  be  realized  in  the  elec¬ 
tronic  systems  that  have  become  so  ubiquitous  in  telecommunications,  computing, 
processing,  and  information  technology.  This  thesis  studies  self-consistent  calcu¬ 
lations  to  understand  how  vacancies  in  the  graphene  monolayer  channel  can  hinder 
such  performance  of  a  predicted  room-temperature  exciton  condensate.  Vacan¬ 
cies  within  a  coherence  length  of  the  contacts  significantly  obstruct  performance 
by  effectively  reducing  the  width  over  which  interlayer  transport  occurs.  For  a 
selected  value  of  4%  vacancies  in  this  region,  tunneling  current  at  a  selected  bias 
can  drop  by  more  than  80%  compared  to  the  ideal  scenario.  The  reduced  width 
of  the  system  caused  by  the  presence  of  the  top  layer  vacancies  produces  a  lin¬ 
ear  dependence  on  the  critical  current  as  the  vacancy  concentration  is  increased. 
Vacancies  outside  of  the  coherence  length  have  little  effect  on  the  interlayer  con¬ 
ductivity,  showing  a  critical  current  with  an  approximate  square  root  dependence 
on  vacancy  concentration  in  the  top  layer.  Critical  current  degrades  up  to  30% 
due  to  a  phase  transition  at  smaller  bias.  The  reduction  is  due  to  scattering  of 
layer  quasiparticles  from  localized  mid-gap  states,  which  modifies  the  average 
interlayer  phase  relationship  between  the  two  layers.  Concentrations  of  vacan¬ 
cies  larger  than  4%  in  one  of  the  layers  prevent  the  condensate  from  forming  in  a 
steady  state. 

The  thesis  thus  provides  a  concise  model  to  show  how  channel  disorder  affects 
superfluid  transport.  As  the  ideal  system  is  difficult  to  realize,  it  provides  exper¬ 
imentalists  a  numerical  framework  for  characterizing  fabricated  systems  domi¬ 
nated  by  short-range  scattering  impurities.  Future  work  includes  adding  the  effect 
of  phonon  scattering  to  the  lattice  to  see  how  superfluid  density  and  critical  tunnel¬ 
ing  current  decay  as  temperature  increases.  Aside  from  lattice  disorder,  phonon 
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coupling  to  the  substrate  is  a  dominant  scattering  mechanism  that  must  be  studied 
in  order  to  have  a  complete  picture  of  room-temperature  performance  degradation 
from  an  ideal  system. 
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